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Comment on "Basis Restricted Transformations 
and Performance Measure for Spectral Representations". 

V. R. ALGAZI 
B. J. FIND 


ABSTRACT : In this correspondence, we comment on some of the shortcoming 

for source encoding problems, of the measure of distance between ortho- 
normal representations proposed by Pearl [1]. 

We further examine, using rate distortion theory, whether rates 
achievable for a given distortion level by the use of basis-restricted 
transformations have some inherent dissimilarity independently of input 
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Comment on "Basis Restricted Transformations and 
Performance Measure for Spectral Representations". 

Pearl [1] has proposed a distance measure to characterize "inherent dissim- 
ilarities" between orthonormal transformations. With this measure the distance 
between orthonormal transformations is "independent of input statistics". One 
of the applications mentioned for basis restricted transformations is to the 
source encoding of random processes. In this correspondence we comment on some 
of the shortcomings of Pearl's distance measure when applied to source encoding 
problems. We further examine, using rate distortion theory, whether rates 
achievable for a given distortion level by the use of basis-restricted trans- 
formations have some inherent dissimilarity independently of input statistics 
and of other implied parameters. 

We note first that dissimilarity betv/een two transformations is defined by 
Pearl for the class of random processes which have a diagonal correlation matrix 
for one of the transformations. The transformation is thus the Karhunen-Loeve 
expansion of a random process adapted to the specific transformation. Under 
these circumstances independence on input statistics achieved by averaging 
over a class of processes cannot be dissociated of the class of processes con- 
sidered. Moreover this construction does not provide any insight into a common 
question of interest in the use of basis-restricted transformations for source 
encoding: Given a random process with known correlation matrix how do different 

transformations rank in the encoding of this random process? This is the specif i 
question that Pearl refers to when he comments on the asymptotic distance between 
the Walsh-Hadamard and the Fourier transformation for large block size. However 
Pearl's distance measure also indicates a diverging distance between the Karhunen 
Loeve (which is then the Kronecker or natural transformation) and Fourier trans- 
formation for large block size. On this basis one cannot feel confident that 
Pearl's distance between Walsh and Fourier transforms provides an answer to a 


comiron conjecture. 


The evaluation of the relative performance of orthonormal transformations 
in source encoding problems can be done by assuming that the expansion coeffi- 
cients are encoded independently and the rate distortion bound evaluated for 
each coefficient. This approach was used by Goblick and Hoi singer [2] to eval- 
uate the performance of the Fourier expansion for source encoding and has been 
extended to other "basis-restricted transformations" by Pearl. 

For example we show in Figure 1 such rate distortion functions for a 
Gaussian Markov random process with correlation function R(n) = (1) 

n = 0, + 1, + 2 --. These curves are for a = + 0.05 and clock size N = 32. 
Similar curves are given by Pearl , Andrev/s and Pratt for N = 256 [3]. Note 
that the Haar and Walsh-Hadamard transformations behave similarly, in contra- 
diction to Pearl's distance. The similarity of performance can be explained 
by exploiting the formal relations between the Haar and Walsh Hadamard trans- 
formations [4]. Note also that for the semi -logarithmic scale used for the 
mean square distortion the rate-distortion per sample curves are parallel 
straight lines. This occurs for mean-square distortions lower than the smallest 
variance of expansion coefficients since we have then 


N 2 

n o] 

R(e) = ^ log2 


N 


( 2 ) 


1 

Therefore dj^^ j ^2 R](e) - ^2^^^ " N ^°^2 ^ ^'^li'^’^2i^ which subscripts 1 

and 2 denote two distinct orthonormal transformations. Mote that dj^^ j ^2 
no longer dependent on the specific e . 

To illustrate the dependence of source statistics we show in figure (2) 
dj^ as a function of a , defined in (1) for the fixed block size N = 32. It 
is noteworthy that the curves for Walsh Hadamard and Fourier transformations 
intersect. Thus no general conclusion can be obtained on the relative merit 



of these two orthonormal transforniations in source encoding problems. Averaging 
on some set of statistics is possible but the use of the results is dependent 
on the specific application. As shown in Figure (3), when the block size in- 
creases, distance d^^ between the Fourier and Karhunen Loeve transforms 
decreases to zero as expected. It can be shown that the distance between the 
Haar and Karhunen-Loeve transforms increases to a limit and it is conjectured 
that the Walsh Hadamard transform has a similar limiting behavior. 
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A UNIFIED TREATMENT OF DISCRETE UNITARY TRANSFORMS 
WITH A FAST ALGORITHM 

by 

i* ■f*i‘ 

Bernard J. Fino and V. Ralph Algazi 


Abstract 

A set of recursive rules which generate unitary transforms with a 
fast algorithm are presented. For each rule, simple relations gives the 
number of elementary operations required by the fast algorithm. The 
common Fourier, Walsh-Hadamard (W-H) , Slant and Haar transforms are 
expressed with these rules. The framework developed allows the intro- 
duction of generalized transforms which include all connnon transforms in 
a large class of "Identical computation transforms." A systematic and 
unified view is provided for unitary transforms which have appeared in the 
literature and for a number of new transforms of potential interest. 
Generalization to complex and multidimensional unitary transforms is 
considered. 
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Introduction 


The dissemination of the Fast Fourier Transform algorithms, 
originally introduced by Good [1] , and knovm as Cooley-Tukey [2] and 
Sande-Tukey [3] algorithms, has resulted in a large extension in the 
range of applications of the well known Fourier transform. Recently the 
Walsh-Hadamard transform, also with a fast algorithm [4] has drawn 
considerable interest [5]. The Haar transform although closely related 
to the Walsh-Hadamard transform [6] and potentially of interest [7], has 
received much less attention. These transforms have been used successfully 
for error free signal representation [8], pattern classification [4], [9], 
speech signal encoding [10] and above all for picture encoding [11] , [12], 
[13], Only a few transforms have been considered in these applications 
while many other transforms could be of interest. Some workers have 
considered the definition of generalized transforms and we mention the 
works by Andrews, et al [14], [15], [16], Rao, et al [17], [18] and 
Harmuth [19]. 

In this paper we present a unified view of discrete unitary trans- 
forms with a fast algorithm. A discrete unitary transform is characterized 
by a unitary matrix [T] such that [T] [T**"] = [I] where * denote conjugate 
"t" transpose and [I] is the unity matrix of same order as [T] , say N. 

For mathematicians a unitary matrix expresses a rotation of the orthonormal 
basis and preserves the Euclidian norm UvU of a vector V. In this N 
dimensional space, Ivl - V.V*^ representation, this property 

means energy conservation and an easy expression of the mean square error 
when some components of the signal are ignored in the new base. The 

computation of the transformed vector W of V by the transform [T] such 

2 

that W = [T]V usually requires N multiplications and N(H-l) additions. 
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For some specific transforms of Interest such as the Fourier, Walsh- 
Hadamard transforms a fast algorithm has been found which requires fewer 
elementary operations. The analysis of these fast algorithms has been 
done by factorization of the matrix [T] into a set of largely sparse 
matrices, each expressing a stage of computation. This is the approach 
followed by Good [1] in his original paper which lead to the Fast Fourier 
Transform [2], [3] the Fast Walsh Transform [4] and other known fast 
transforms. 

Here we consider recursive rules^ for the generation of unitary 
transforms having a fast algorithm. These rules allow us to generate 
large classes of such transforms, many of which are new and possibly of 
interest, and to give general formulas for the number of elementary 
operations required by the corresponding fast algorithm. 

1. Recursive Generative Rules : 

We shall present three rules which generate a new unitary matrix 
from some original unitary matrices. For each rule we relate the number 
of elementary operations for the new transform to the number of elementary 
operations of the same type required by the original transforms. For rule 1 
there is only one original matrix, for rule 2 two, and for rule 3 a set of 
original matrices. 


we denote by "rule" a set of operations performed in a prescribed order. 
We reserve the term "operation" for the elementary operations such as 
additions, multiplications, etc. which determine the computational 
complexity of a transform. 
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Rule 1; Operations on the columns of a unitary matrix ; 

Given a unitary matrix [T] , two obvious operations on the columns 
yield another unitary matrix of some order: 

a) permutation of the columns: This operation does not require 

any computation. In the computational process, this operation can be 
performed by applying the permutation to the coefficients of the input 
vector Instead of the columns themselves. 

b) multiplication of a column by a root of unity: This operation 

introduces a complex multiplication if the root of unity is not + 1 or 

+ j (j = (see footnote 2). 

These operations on the columns may be expressed by a matrix 

i 0 4 

product [T] [D] with [D] such that = e“^ ^ if column k is to be replaced 
by column 1 multiplied by the root of unity, e and all other entries 

of [D] are null. 

Rule 2 ; Rotation of rows by a unitary matrix 

Consider a unitary matrix [T] of order N. The N row vectors form 

an orthonormal basis for the N dimensional space they span, m rows 

vectors of [T] form an orthonormal basis for subspace S^. If these m 

vectors are rotated by a unitary matrix [U] of order m, we obtain a new 

orthonormal basis^^ for S . The remaining unchanged N-m rows of [T] are 

m 

an orthonormal basis of the subspace orthogonal to and form with 

a new orthonormal basis for . Thus, the matrix [T'] obtained after 
rotation of the m rows by the unitary matrix [U] is unitary. 

2 

Multiplications by +1 ^^d ^ may be counted as operations if the hardware 
realization of the algorithm is not able to keep track of them. However, 
for the error analysis of the algorithm these multiplications, even if 
they are performed, do not introduce any error. 
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Some particular cases of interest are: 

a) multiplication of the whole matrix by a unitary matrix of 
the same order 

b) permutation of the rows (multiplication by a peirmutation 
matrix) 

c) multiplication of a row by any root of unity. 

The operation b) and c) can be represented by the matrix product 
[D] [T] where [D] is, as before, such that = e^^^ if row i of T is 

replaced by row k multiplied by the root of unity, e“^ and all other 
entries of [D] are null. 


Number of Elementary Operations; 

If transforms T and U require respectively t and u elementary 
operations of a specific kind, it is obvious that the transform T’ will 
require at most t' of these operations with 


t' = t + u (1) 

(It may happen that [T*] so generated has a simpler algorithm). 

Equation (1) applies Independently to any type of elementary 
operation, additions, real and complex multiplications as well as any 
other specific operation (e.g. shift, multiplication by /J . . . etc.) 


Rule 3: Generalized Kronecker Product: 

Given two sets of unitary matrices, set of m matrices [ j 
(i=0, . . . ,m-l) all of order n and set I" of n matrices | B^j (i=0 , . . . ,n-l) , 
all of order m, we define the generalized Kronecker product of the sets 


and {CB}. denoted 


to be the square matrix [C] of order 


(mn) such that 
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( 2 ) 


w u ' 

C. . = C . » = A I • B , 

1 , 2 um+v , u urn? uu ww 


with i = um+w 


u. 

u'=0 , ...» n-1 

j = u’mhw' 


w. 

w'=0 , . . . , m-1 

a) [C] is a unitary 

matrix: 



Proof : 




mn-1 * 

n-1 

m-1 


2 ^jk 

' E 

E 

c c 

um+w , vmhz u ’ mfw ' , vmf z 

k=0 

y=0 

o 

II 

N 


with k = vm +2 

v=0 

9 * * • 9 

n-1 


z=0 , . . . , m-1 


Using (2) 


mn-1 


E 

k=0 


^ik ^jk 


n-1 

m-1 

a'^ 

uv 

4* .. 

= E 

E 

. * w 
A r 

U V 

v=0 

N 

II 

o 



n-1 

A ” 
uv 


m-1 

-E 

A* W 

A t 

U V 

E 

v=0 



z=0 


„v „*v 


•5^ by orthonormality of 


w 


n-1 


= 6 


w,w' ^ 
v=0 


a" a* 

uv u V 



by orthonormality of 


= <5 ,6 , = 6 , . 

ww uu ij 


where 6^^ is the Kronecker delta 6^^ = 1 if i=j 


= 0 otherwise. 

This proves that [C] is a unitary matrix. 

In the particular case in which the matrices j^A^j “ ^aJ 


QED 

are all 

I J L J 

identical, and also the matrices then the generalized Kronecker 
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product 0 reduces to the usual Kronecker product of matrices 

[14]: [A]0[B]. 

b) Factorization of [C] i We now prove that 

[C] = [P^] [Dlag{^}] [P] [Dlag{Cg}] (3) 

where [DlagO^)] and [Dlag{C0}] are block diagonal matrices formed with 


the matrices of the sets } and 

permutation matrix of order mn such that 

S, = v'm+z' and z' , v=0, m-1; z, v* = 

a generalization of the factorization of 
Good matrices [14] . 

(see Fig. 1) and [P] is the 

P, „ = 6 ,6 . when k = vn+z, 

ki, vz' zv' 

=0, n-1. Equation (3) is 

a simple Kronecker product into 

Proof : 


[Diag04}]^,j^ = 

with k' = v"n+z" 

[Diag®]^^ - 6^,^, B^,^, 

j = u*m + w’ 

^ik' " ‘^uz" ‘^wv" 

i = um + w 


We evaluate an element of the matrix on the. right hand side of (3) 


mn-1 

mn-1 

mn-1 


s 

z 

z 

[p''] [DiagU}]j^,j^ = 

o 

II 

k=0 

1=0 


m-1 n-1 

m-1 

n-1 

n-1 m-1 


Z 

d d „ d „ a'^,, ,S ,6 , 6 , , b“, 

y ^ /—J! uz wv w z z vz zv’ u'v' Z w 

<_ 

if 

o 

N 

if 

O 

< 

II 

o 

o 

II 

N 

<! 

II 

o 

N 

11 

o 

^ , 

uu 

b“', 

ww 

-=1J 

QED. 
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c) Number of elementary operations ; 

With the computational blocks corresponding to the transforms 
A°, ^ and b” the factorization of equation (3) 

leads directly to the computational block of the transform C given in 
Figure 2 . 

From the structure of the algorithm of Figure 2 it is easy to see 
that if the matrices [A^] (i=0, •••, m-1) and [b'^] (j=0 , . . . ,n-l) 
have algorithms requiring respectively p^ and q^ elementary 
operations of a specific type, their generalized Kronecker product {CJ 
will require of these operations with 


m-1 n-1 



i=0 j=0 


In the particular case of a simple Kronecker product p^ = p^ and 


P = m P + n q (5) 

mn n ^m 


Note that the use of rule 1 and 2 only increases the number of elementary 

operations while the order of the generated transform does not change. 

For rule 3, even if [A^] and [B^] do not have fast algorithms and thus 
2 2 

require n and m elementary operations, [C] requires a maximum of 
2 

(m+n)mn ^ (mn) (for m, n > 1) elementary operations. 

The results of equations (1), (4) and (5) are important: for every 
transform generated with the recursive rules presented, they give a 
simple and systematic way to estimate its computational complexity. 
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2. Identical Computation (IC) Family 
The generative rules defined above create a unified framework 
for the known fast unitary transforms, introduce new transforms, and 
allow assessment of the computational complexity of such transforms. 
In this paper, one large family of transforms is considered; the 


"identical computation transforms" that we discuss now. 

We denote by {^}@ generalized Kronecker product of 

a set of q matrices [A^] (K=0 , ..., q-1) of order p and a set 

P 

of p identical matrices [B^] of order q. [B^] will be called a core 
\c 

matrix and [A ] a parent matrix. The IC transforms are recursively 
P 

generated from a unique class, (3> parent matrices of some order f 
and an original core matrix [©] of order q. An IC transform of order 
(qf"^) is then obtained from the original core matrix [©] by the recursive 


formulas ; 


(IC^fl - IDqjl [{j4) ® [OJ] [D'j] 

[IC 1 - [D [IC [®' J 

qf qf qf qf 


( 6 ) 


where the matrices [D] and [D'] express respectively a reordering 
followed by multiplications by roots, of unity of the rows and the columns. 
All parent matrices of .... belong to C. 

The common characteristic of all the transforms of the IC family is 
that their algorithms only use in any computation intermediate results 
obtained from the input vector through identical computations (so the name 
of the family). This property provides a uniform treatment of successive 
components of the input vector if we consider that any parent matrix 
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treats uniformly its input vector. For this family, all the normalizations 
can be delayed to the last stage of computation. 

We shall consider different choices for the original matrix [©], 
the class Q of parent matrices, the matrices [D] and [D'] and the sets 
We first show that the basic transforms, Fourier* Walsh-Had amard 
and Haar, are IC transforms. 

3. Basic Transforms; Fourier, Walsh-Hadamard , Haar ; 

In this section with the help of the generative rules, we examine 
the well known Fourier, W-H , and Haar transforms. This approach allows 
the derivation of some new results concerning the number of multipli- 
cations required by a FFT of composite order, a concise presentation of 
the different definitions of the W-H transform, and simple definitions 
of the Haar transform. In addition it makes apparent the common 
structure of these transforms. This will lead in the next section to 
the definition of families of transforms between the basic transforms. 

In the following we emphasize specific orderings for the basic 

3 

transforms: frequencies for the Fourier transform, zequencies for the 
W-H transform and rank for the Haar transform. These orderings have 
proved to be useful in signal encoding because they concentrate the 
signal energy into the first transform coefficients, for some image 
models [21] . 


This terminology has been introduced by Yuen I20J. The zequency is the 
number of zero crossings. 
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3.1 Generalized Fast Fourier Transform of Composite Order 
a) decomposition theorem ; 

Given the Fourier matrices [F ] and [F ] of orders p 

P q 

respectively, the matrix [^pq] such that 
[Fpq] = [{[F^]}® [Fp]] [P]’" 


is the Fourier matrix of order pq. The set {[F ]}, k=0, ... p-1 

q 

such that 


[F ] = [F ] [D ] where 
q q k 


Vn * 

[Dj^] is a diagonal matrix such that (Dj^) = exp(-2irj ) 


.. 


u'u' 


[P] is the permutation matrix such that 


pq 


P = 6 6, with s 

St uz kw 


= uq + k t - wp + z 


и, z < p 

к, w < q 


Proof; We denote [{[F^]}® ^^pq^ 


^^pq^ ug + k, u'g + k' ^^p^uu'* ^^q^kk' 


-2tt3 


= (F ) . e 

uu 


ku’ 

pq 




„ , ,uu ' , ku ' , kk \ 

T -2ttj ( + + ) 

L- e p pq q 


[F„J = [F* ] • [P]^ (F ) = (F’ ) .6 

P^ ug+k, wp+z uq+k, u*q+k- 


and q 


(7) 

is 


C8) 


' "^wk’ 
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= e 


_2Trj(H5.+ M. + ^) _2Trj(lia±kLJ5^2±El 

p pq q _ „ pq 


QED. 


= e 


Note that I^pq3 symmetric and orthonormal so that 
—1 * 

([Fpq] ) = u.se of (3) and (7) it is possible to derive a 

new expression for [F ] : 

pq 


[Fpq] = [[Fp] 0 {[F^*"]}] [P] 


(9) 


with [F^^] = [Dj^] fF^] 


If [Fp] and [F^] require respectively c_yl:p and complex additions, 
(jtip and^^ complex multiplications, will require by application 
of (4). 


where C 

p,q 


PkAq + q 'Jp 

complex additions 

(10) 

p;_y‘l/t + qU^ + C 

q p p>q 

complex multiplications 

(11) 


is the number of complex multiplications introduced by (8) . 


C 

Pfq 


c 


p,q 


pq if all the factors including + + j considered. 

(p-1) (q-1) if the factors + 1 are discarded 
(p-1) (q-l)“l if the factors + j are also discarded 
(when (pq) is a power of 2). 


b) Generalized FFT of composite order 

If the order of the Fourier transform is composite, i.e. N = . . . 

p^, the previous decomposition theorem yields the well known FFT 


- 11 - 



algorithms [2] [3] detailed by Classman [22] in the most general case. 
The recursive use of the formulae (10) and (11) gives the number of 
required operations. In the case of N = r^ we can solve these recursive 
equations: this is the case of FFT of radix r. 


oA = r’^ ^ + r cA or = nr” ^ 

r” ' ^ r” 

( 12 ) 

iJ[j[ = r” ^u\JL + tlAA + (r-ct) (r” ^-ct)- 6 or 
n r n— i 

r r 


n-1 


= n r” ^'jU^ + (r-a) [(n-1) r” ^ - a ] - 


”"^-1 

(r 1) 

r-1 


(13) 


(a, g depend on the value of C ) 

p,q 

The radices 2, 4, 8 and 16 have been considered in the literature. 

For the radix 2, which gives the most popular FFT, the recursive 
relations given by the decomposition theorem are 


[F^^] = [{[F 2 ]}® [F^n-1^^ 


(14) 


and [F ] = [[F {[F [P] 

2 ” 2”’^ ^ 


wa 


ith [p 


ki 1 




“] 

1 e 

, 'k, 1 

1 1 ~ 

-2lTj — 

1 2 ” 

1 -e 

and [F 2 ] = ^ 

-2irj — -2nj — 
e 2 " ® 2 ” 


(15) 


The algorithm corresponding to the recursive formula (14) and 
obtained by recursive use of Fig. 2 is shown in Fig. 3a; it can be 
arranged equivalently with all operations "in place" as shown in Fig. 3b 
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which is the classical diagram of the Cooley-Tuckey [2] algorithm 
with decimation in time. 

The algorithm corresponding to the formula (15) is the Sande- 
Tukey [3] algorithm with decimation in frequency and is shown in Fig. 3- 
c and d. 

For these Figures the factors are 

a = 1 
o 

a^ = exp(-27Tj/8) 
a^ = exp(-47Tj/8) = -j 
a^ = exp(-6-trj /8) 

We can compare the FFT with radices 2,4,8 and 16 for transforms 
n 

n 

of order N « 2 = r (n is then a multiple of 12) . The formulas 

(12) and (13) then give: 



2 

The column_^,|^| has been given by Singleton [23] . In fact our approach 
r 

allows the evaluation of the computational complexity for any composite 
order, in particular for mixed radix FFT. 
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The factors + 1 are easy to track in the algorithms and for most 

realizations multiplications by + 1 are not performed. The factors 

+ j appear in various places in the algorithms and in most realizations 

multiplications by + j are performed; however, in an error analysis 

these multiplications do not introduce any rounding error and the 

column ylyP is then of interest, 
n 

r • 

3.2 Walsh-Hadamard Transform ; 

The W-H functions are well known and the results presented in this 
section are explicit or implicit in many publications. Here, we 
wish to express these results in terms of our generative rules; we 
think that the following compact notation clarifies the relations 
between the various orderings of the W-H functions and the different 
algorithms. This approach also makes apparent the common structure of 
the W-H transform with the Fourier and Haar transforms. 

T^ee distinct orderings or rows of the W-H transform are commonly 
used and are of interest (see the discussion by Yuen [20]). For each of 
these orderings there exists a recursive matrix definition: 

a) " natural order" It is obtained by simple Kronecker product 
without any permutation 

[WH nat.] = [F [WH nat.] (16) 

with the original core matrix This relation gives directly, by 

recursive use of Fig. 2, the fast algorithm of Fig. 4a (without the 
reorderings) . A different presentation of this algorithm with identical 
stages of computation is given in Fig. 4b. 
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b) Paley's ordering ; Used originally by Paley {24] it seems 
more suitable for mathematical developments than the other orderings. 

The recursive relations introduced by Yuen [20] are expressed by the 
matrix relation 

[WH pal.] = [[F_]0 [WH ^ ^ pal.]] [P]^ (17) 

2 2 

with the original core matrix [F 2 ] and [P] as defined previously for 
the Fourier transform (see section 3.1). This relation gives the 
algorithm of Fig. 4c (without the bit-inversion reordering) by recursive 
use of Fig. 2. A different presentation of this algorithm can be 
found in [6] . 

c) zequency ordering ; This is the original ordering by Walsh 
[25] and is the ordering of interest for signal encoding because it 
ranks the transform coefficients roughly according to their variances 
for signal statistics commonly encountered in practice. The generating 
process*^-bU > introduced in [6], defines recursively the W-H matrices in 
zequency order. We can express it by the matrix relation 

zeq.] = [W] [[F 2 ]® 

where [W] denotes the reordering of the process A zequency ordered 

algorithm has been investigated by Manz [26]. 

Although the zequency ordering could be generated recursively, the 
corresponding algorithm would not be simple and it may be preferable to 
obtain the W-H transform in natural or Paley's order and then perform a 
global reordering. 
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Given the transform coefficients in Paley's order a 
"bit inversion" reordering, denoted in matrix form by [Q] , is necessary 
to put them in zequency order: in a bit-inversion permutation, 

consecutive coefficients with bit of the binary representation 
of their indexes equal to 1 are put in reverse order. This operation 
is performed for all bits starting from the rightmost bit. With 8 
coefficients to reorder this means the following permutation: 


index of 
coefficients in 
Paley's order 

binary 

representation 

reverse order 
for middle 
bit 

reverse order 
for leftmost 
bit 

final 

order 

(see Fig. 4c) 

0 

000 

^ 000 

000 

0 

1 

001 

^ 001 

► 001 

1 

2 

010 

Oil 

— ^ on 

3 

3 

Oil 

010 

► 010 

2 

4 

100 

► 100 

, no 

6 

5 

101 

— 101 \ 

A 

7 

6 

110 .. 

^ 111 

101 

5 

7 

111 

no 

A 

100 

4 


Given the transform coefficients in natural order, a bit reversal 
ordering, denoted by the matrix [R] , will order them in Paley's order and 
a bit-inversion will order them in zequency order. For 8 coefficients 
to reorder we have: 
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in any of the three orderings, these reorderings can be performed on 
the columns as well as on the rows. Hence we have the matrix relations; 


[WH zeq.] = [Q] [WH pal.] = [T-JH pal.] [Q]'" 

2" 2" 2 


[WH zeq.] = [Q] [R] [WH nat.] = [WH nat.J [R]'^ [Q]*^ 
2 2 2 


Since the W-H matrices are their own inverses, we have also using Cl6) 
and (17) the following recursive relations 


[WH nat.] = [P]*^ [[WH , nat.]® [F ]] [P] 
2 ^^ 2 ^^ ^ ^ 


[WH ^ pal.] = [[WH pal.]@ [F 2 ]] [P] 


These relations however do not give different algorithms. All these 
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algorithms differ only by reordering and so have the same number of 
additions given by 


■An-^ • An-l'"'- 'A ' ^ 


which gives 




= n2“, a well known result. 

- 2” 


3.3 Haar transform ; 

The Haar transform is usually defined from the Haar functions 111] 


The Haar matrix 

of 

order 8 

[Hg] 

ordered 

by 

ranks is 

as 


1 

1 

1 

1 

1 

1 

1 

1 


1 

1 

1 

1 

-1 

-1 

-1 

-1 






0 

0 

0 

0 

1 

0 

0 

0 

0 

x/T 


-yr 



2 

-2 

0 

0 

0 

0 

0 

0 


0 

0 

2 

-2 

0 

0 

0 

0 


0 

0 

0 

0 

2 

-2 

0 

0 


0 

0 

0 

0 

0 

0 

2 

-2 




Zones 

0 




Here we use the generative rules to define recursively the Haar matrices 
and we have found two definitions : 

1) The Haar matrix of order 2^^ is obtained from the Haar matrix of 

1 0 


order 2^^ ^ by simple Kronecker product with. JI^J =. 


0 1 


followed by 


rotation of the rows 0 and 2^ by is the process of [6J , 
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The rows are obtained in "natural” order. To reorder them 

by their ranks, we need a "zonal bit reversal" ordering. A 

zone as defined in [6] is a set of coefficients with indexes between two 

successive powers of 2. A "zonal bit reversal" ordering is a bit-reversal 

followed by a reordering in the original order inside each zone. For 8 

coefficients the zonal bit reversal ordering gives; 

Index Binary representation Bit reversal Reordering final order 

inside zones (see Fig. 5a) 

0 000 — 000 } 000 0 



With both definitions we obtain by recursive application of the diagram 
of Fig. 2 the algorithm of Fig. 5a. This algorithm can be more conveniently 
organized as shown in Fig. 5b and give the rows directly ordered by their 
rank. 

By application of (4) we obtain the following recursive formula for 
the number of additions: 
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J{ = 2 ► f + 2. Hence = 2(2"-l) with /jl. = 2. 2^ ^ 

2" 2^~^ 2*^ ^ 

normalizations are also required. A modified Haar transform obtained from 
the Haar transform by permutation of its columns is related to the 
Fourier transform (see later section 7.3): it is defined recursively by; 

[MH^^] = [Z] [{[F2], [I2]. .... [I2]}® [P]^ (20) 


Globally the permutations [Zj perform a bit-reversal ordering inside each 
zone . 

The modified Haar matrix of order 8 is as follows 


'“>81 yr 


-1 


\ji 0 

0 72 


0 0 


-1 


-1 


0 0 


0 -2 


-2 


0 0 0 -2 


-1 


-72" 0 0 -yi 0 

Q -yfi 0 0 -yr 


0 0 


and its algorithm is given in Fig. 5c. 


4. Generalizations of the basic transforms 

For the three basic transforms we have found recursive definitions 
with a matrix formula similar to (6) . The direct comparison of these 
definitions suggests the generalization of the basic transforms to families 
between them. The simplest generalizations are two families between the 
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Fourier and W-H transforms, and between the W-H and Haar transforms. A 
larger generalization is the IC 2 family which includes all three basic 
transforms with parent matrices of order 2. A number of these generalized 
transforms has recently been discussed independently. We would like to 
show that they fall easily within the framework we have developed and that 
further generalizations are clearly possible. 


4.1. Family between W-H. and F ; 

If we compare the recursive generation of the Fourier matrices with 
radix 2 (14) and the W-H transform (17) , we notice that they differ only 
by a set of factors. If V7e exclude the reordering of the rows, we see 
that we can easily generalize the W-H and Fourier transforms to a large 
family of unitary transforms given by the recursive formula 






IPJ 


where [F2(9)] 



1 

1 


exp(-j0) 
-exp(-j0 ) 


This family includes the W-H and Fourier transforms for the appropriate 
choices of the parameters 0^ . . . 0 ' . 

Two families have appeared in the literature for special choices of 
these parameters: 

1) 0, = v/here c is a real 

scalar varying from 0 (W-H transform) to 1 (Fourier transform) . The 
corresponding transform has been called "general spectral analyzer" I14J, 
[15]. 
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2it k 


2) e^. 


if k mod (2^^ = 0 


where g is an integer varying from 0 (W-H transform) to n-1 (Fourier 

transform), and 0^^ = 0 otherwise. The corresponding transform is the 

"Generalized Discrete Transform" [17]; we will denote it [GT®^] . 

2 ^ 

Fig. 6 shows the matrix [GT ] and its fast algorithm. Many other choices 

Id 

for these factors are obviously possible and these two special choices 
do not seem to bear any exceptional importance. 

As an example of use of our general formulas ve compute now the 
required number of multiplications for the Generalized 
Discrete Transform. There are 2^-1 factors different from 1. 


So j\j^ = 2 j\j[^ T + 2^-1 for n - 1 > g 
„n .,n-l 


otherwise. 


So that = g 2^ ^ - 2® + 1 


If we do not count the multiplications by jj , we find similarly 

u\.[\ = (8-1) 2""^-2S + 2 
2 

4.2 Family between Hear and W-H 

In [6] one of the authors has presented a family of transforms between 
the Haar and W-H transforms. This family was obtained by replacing W-H 
transforms of lower order by Haar transforms in the decomposition of the 
fast algorithm of a W-H transform. Now we have further decomposed the fast 
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algorithms of the two transforms up to similar recursive formulas (16) 
and (19) or (17) and (20). Obviously, if we choose any of the 2^ ^ 
parent matrices needed to generate the matrix of order 2^ to be either 
[F 2 ] or [I 2 ] I we obtain a large family of unitary matrices which includes 
the Haar, W-H, and unity matrices. There are 




+ 2 


n-1 


= 2 


2’^-l 


members of order 2^ in this family. 

The number of additions is obviously twice the number of parent 
matrices equal to [F 2 ]. For the W-H transform we have n2^ ^ such matrices 
and therefore n2^ additions. For the Haar transform we have 2^^ ^ + 2^ ^ 

... + 1 = 2^-1 such matrices [F 2 ] and so 2(2’^-!) additions. The number of 
normalizations varies from O(W-H) to 2’^-2 (the normalizing factors come by 
pairs and all pairs are different in the v/orst case) . No multiplication is 
required during the computation. At the order 8, 2^ matrices are in the 
family. We show in Figure 7 one of them with its fast algorithm. 

Assume as a particular case that we choose the parent matrices of 
the recursive formulas to be 

[F^] k = 0, .... 2^"^-! 

with 

[F^] = [F 2 ]for k = 0 mod (2P'''^“") 

Ic 

F = [I 2 ] otherwise 


where p is the stage of computation up to n when we generate a transform of 
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order 2 ^ and h an index lower than n. • 

Then, i-f the recursive formula used is similar to (16) and (19) we 
obtain a subclass of n transforms; for h = 0 we have the W-H transform and 
for h =< n - 1 we have the Haar transform, both in natural order. 

If the recursive formula is similar to (17) and (20) (with the per- 
mutation [P]*" of the columns) we still obtain n transforms: for h = 0 we 

have the W-H transform in Paley’s order and for h = n - 1 we have an 

It 

unordered modified Haar transform. We denote these transforms [WHH^n]. 

2 ^ 

Fig. 8 shows [WHHj^g] and its fast algorithm. 

4.3 IC ^ family ; 

To generate the family between W-H and W we have introduced a set of 
factors into the recursive formula for the W-H transform. To generate the 
family between W-H and Haar, we have replaced some parent matrices by the 
identity matrix [I 2 J in the same recursive definition of the W-H transform. 
If we allow simultaneously both operations we generate a larger family 
that we call IC 2 • 

More formally if [T , 1 is a member of IC^ of order 2^^ a member 
of order 2^^ is given by 


[T^„] = [D^] [{[CqJ, ..., 


fCn_i 

2 -1 2 


n-1 


JJ 


where [D^] and [^ 2 ] are permutation matrices and C^, 



are either 


[I,] or [F re)] = i 
^ ^ /2 


1 

1 


exp (-he ) 
-exp (-je ) 


- 24 - 




(4) This notion of degree of freedom is an extension of a concept introduced 
by Andrews and Caspari [16]. For them the degree of freedom of a class of 
matrices is the number of free parameters required to define this class. 

This definition is ambiguous when the constraints which define a class cannot 
be reduce^si -.t o a set of free parameters. For example, the unitary matrices 
of order 2 are given 1 degree of freedom in [16] when in fact, on the real 
numbers, the most general matrix is 

cos a sin a e = + 1 

e sin a -e cos a a ^ [0,2tt] 

and on the complex numbers the general solution depends on 4 angles 
G [0,2it] and 2 binary choices. Our approach is to track as far as possible 
the reduction to independent choices. If it can be reduced to a number of 
free parameters our degree of freedom will be the number of these parameters. 
Note that the relations (1) and (4) apply also to the recursive computation 
of the degree of freedom of a class. Note also that the degree of freedom 
has generally no relation with the computational complexity (which varies 
usually for the transforms of a class) . 
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where is the number of parent matrices different from [I^] at this 
stage, and the number of factors different from +1 (and maybe + j ) at 
this stage. 

For Haar A, = 1 for any k and = 2(2”-l) 


For W-H 


= and = n2^ 


For Fourier = n2”^ and L, = 2^ 2^ ^-1 or 2^ ^-2, which yield 

2^ k 


,k-2 


the results of section 2.3 (radix 2). 

We present now an example of interest in the IC 2 family; a class of 
transforms which make a discrete transition between the 3 basic transforms 
and which we call therefore the WFH class. 

Each transform of this class is indexed by two positive integer para- 
meters h and g such that h + g < n when 2^^ is the order of the transform 

2 ll 

and is denoted fWFH®’ ]. 


[WFH 




n 


] is obtained recursively as the Fourier transform of radix 


^ k 1 

2 (formula 1^ ) but with the parent matrices [P-] k = 0, 2^ such that 


[P^J = [F (2uk/2P)] for k = 0 (mod 2^ ® 

[P^l = [^ 2 ^ k = 0 (mod 2^^^ but 

k ^ 0 (mod 2^ ^ 

[P^] =[ 12 ] otherwise 

where p is the level of computation up to n. 
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We can represent then WFH transforms on a (g, h) plane as shown In 
Fig. 9. With appropriate permutation matrices, for h = 0 we have the n 
Generalized Discrete Transforms (see section 4.1), for g = 0 the n WHH 
transforms (see section 4.2). WFH^’^ is the W-H transform, WHH^’ ^ ^ the 
Modified Haar transform and WHh” the Fourier transform. For h + g = 
n - 1 we have a set of n transforms in between the Fourier and Haar trans- 
forms which have been called the Modified Generalized Discrete Transforms and 
defined after much work in [18J. 


5. Other IC transforms ; 

Except for the Fourier transform, we have restricted ourselves 
so far to IC transforms obtained from original core matrix 1^2^ parent 

matrices of order 2. The generative rules have given a unified approach 
of the usual unitary transforms. We now consider some examples with a 
different original core matrix and parent matrices of higher orders. 

The matrices of order 2 are of practical interest for the fast 
algorithm as long as we perform the required operations (specially 
additions) with only two operands at a time. If fast additions involving, 
let us say, f operands, become available, the transforms with parent 
matrices of order f may be of interest. 

Most of the recursive structures of the transforms presented in the 
previous sections can be applied to parent matrices of higher orders 
than 2 . We now give some examples : 

a) different original core matrix 

In the definition of the W-H transform the original core matrix 


[F^J may be replaced by the core matrix 


cos 0 sin 0 


sin 0 -cos 6 


and we obtain a 
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transfom considered by Andrews et al. [15] . This original core matrix 
can be used for all the recursive definitions considered, 
b) Generalized 2 and 3 valued transforms ; 

In the definition of the W-H transform the role of [F 2 ] as original 
core matrix and parent matrix can be performed by any unitary matrix 
[U] of order f. If [U] is an Hadamard matrix (its entries are +1/,^/?) 
the generated matrix of order f’^ will also be an Hadamard matrix. These 
matrices have been called "generalized 2-valued transforms” [19]. 

Similarly we can replace [F 2 ] in the definition of the Haar transform by 
the same matrix [U] and we will generate a unitary matrix with entries 0 
or + 1/C^ where is the normalizing factor of the i*"^ row: these 

matrices are the "generalized 3-valued transforms" [19] . More generally [U] 
can replace (F 2 ] in the definition of the family of transforms between 
W-H and Haar. 


c) IC ^ family : 

The IC 2 family was based on the set the parent matrices. We 

can define similarly the IC^ family based on the class of parent matrices 

of order f which contains [I^] and •••> where column 

of •••» column of [F^] x exp(e ([F^] is the 

f ^-1 

Fourier matrix of order f) . The family IC^ has independent parent 

n f”-i 

matrices chosen in we say that IC^ has y ~ degrees of freedom over 
The required number of additions and multiplication is computed 
recursively as done for (21) and (22) : 




multiplications: xJ{J[ = 




n 


n 


k"l k-2 


f"-\, 


( 24 ) 


where is the number of parent matrices different from [i^] at the k stage 

and Lj^ the number of column multiplications with factors different from 

th 

+ 1 (and mayhe + j) at this k stage. 
d) subfamily of IC^ 

£ 

By analogy to the WFH”’ subfamily of IC- we can define the subfamily 
[WFH®’^] of IC^ as follows: 

fU 

£ tx 

[WFH®* ] is obtained by successive generalized Krbnecker products with 


pP-1. 


0 k f*^ -1 

the sets { [M ] , ... [IT] . . , , [M J) of parent matrices such that 


= [fj] with column i of “ column i of [F^J x 

for k = 0 mod (f^ ^ 


= [F^] for k = 0 mod and 

k ii^ 0 mod (f^ ® 


[M J = [I^J otherwise 


and at each level the permutation matrix [P] , P ^ = 6 6, with s = uf + k, 

St uz kw ’ 

t = wf^ ^ + z, is applied to reorder the columns. 

It is. easy to see that [WFh’^ is the usual Fourier transform of 


order f”: the matrices [WFH^’^I and [WFH^’”'^] have been introduced in the 

_n ..n 
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literature respectively by Chrestenson [27] and V7atari [28], For these 
2 matrices (23) and (24) reduce to the same recursive formulas and denoting 
by ^ the number of additions or multiplications: 


for WFH ’ 
f^ 

for 

,-n 


cp = 


^P„- 


n 


n 




n-1 


f f-1 


(This last result corrects the result given in [15] , page 20J 

The other matrices of the family can be represented in the g-h diagram of 

Fig. 9. 

6. Slant transform : 

The Slant transform has been proposed by Enomoto et al. [29] for the 
order 8. Pratt et al. [30] have generalized this transform to any order 
2*^ and compared its performance with other transforms [31] . In this section 
we want to express the recursive generation of the Slant transform with our 
generative rules and compute the number of elementary operations required 
by its fast algorithm. 

The Slant transforms of orders 4 and 8, [S, ] and [S_] , are as follows 

H O 

(in "natural" order) . 


IM 


A 


1 

-3 

1 

-1 


1 

3 

-1 

-1 


1 

-1 

-3 

1 


X 


1//5 

l/y/5 


Zequencies 

0 

3 

1 

2 
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The rows can be reordered by zequencies with the same permutation as 

the W-H transform in natural order. 

The Slant transform of order 2” in natural order is obtained from the 

n*”l 

Slant transform of order 2 in natural order by simple Kronecker product 
with [F 2 ] followed by rotation of the rows 2^ ^ and 2^^ ^ by the matrix 




- 31 - 



with components linearly decreasing: 


(2 - 1) - 2i 


2^(2^n_ 1) 


But some normalizations can be delayed to the last stage of computation 
n 2 XI X 

and the rows 2 and 2 are rotated by the matrix 


,n-l 


(22n-2_i) 


,n-l 


requiring 2 shifts, 2 additions, 1 multiplication. The corresponding 
algorithm is shown in Fig. 10a. 

Number of elementary operations ; 

Formulas (1) and (5) give; 
for the number of additions: 


^ =2,1 ,+2’^ ^.2 + 2 with = 2 


hence <1 = (n+1) 2^-2 


for the number of shifts: 


S „ - 2 • S^„-l 2 fi2 - » 


hence 2 =2*^-2 

2" 
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for the number of multiplications: 



n-1 


+ 1 


with;_^Ay^ = 0 


hence oA/| = 2^ ^-1 

2 “ 

Finally 2^^ - 2*^ ^-1 normalizations are required at the last stage of 
computation. 

However the algorithm at the order 4 can be performed with 8 
additions, 2 multiplications as shown in Fig. 10b [30] instead of 10 additions 
and 2 shifts as above. The formulas (1) and (5) give then: 
additions : 

e/ = 2 Jk' + 2 "“^ .2 + 2 withtyl! = 8 

2 " ^ 


hence (_A = (2n+l) 2^”^-2 = , 1 _"2 

/.»n 


n-1 


shifts I 

5,1 - 2 • Sl_i + 2 Si - 0 
2 2 

hence 2 = 2^^ ^-2 = 2 “ 2^ ^ 

2 " 2 

multiplications : 

hence M ' = 32’^“^-! =u\j[ + 2^“^ 

2 " 2 ” 

n n— 2 

and as before 2 - 2‘ - 1 normalizations. 
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7. Additional properties and generalizations of unitary transforms ; 

In this section we discuss briefly the complex extension of a real 
transform. We also point out some additional relations between transforms 
suggested by the unified framework presented. 


7.1 Complex extension of a real transform ; 

From a real unitary matrix [RT] with rows RTq» •••» RT^^ -|^> we construct 

a complex extension noted [CT] with rows CTq,...., by creating two 

complex rows CT and CT from two real rows RT RT as follows 
p q m n 


CT = ^ (RT -j RT ) 
p yr m n 


CT = (RT +j RT ) 


(25) 


Then the complex transform^ of a complex input vector 

V = R + j 1 is expressed uniquely from the real transforms of and 
denoted*^ and -5^; 


^p ^n^ ^^'^m "^n^ similarly 


With these relations the properties of complex transforms can be 
deduced from those of the real transform. In the literature, besides the 
real and complex Fourier transforms, the complex W-E transfoms (also 
called Complex BIFORE transform) [32] [33] complex Haar transform (also 
called Complex Modified BIFORE transform) [34] have been defined. 
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Note that the complex W-H transform obtained by relations (25) would have 
entries — ^ — : commonly the rows are then rotated by to give a trans- 

/2 n 

form with entries +1 and jj . The rov/s of the complex W-H transform can 
be ordered according to a generalized frequency defined as the number of 
clockwise rotations around the origin when following cyclically the entries 
of a row. 


7.2 Multidimensional transforms 

The techniques presented for the one dimensional transforms extend 
to multidimensional separable transforms. Let us denote an input array 


of p dimensions by A. , .. .., ^ and the p-dimensional separable transform 

P . 

Then the 

P P 


by T „ . . = T^ . T^ . ... T^ . . 

Uf > . . . . , Up , , . « . . ) ip Uj^ 1^ U 2 ^2 ^ 


transformed array 


x„ X 

12 p 


A . , T . . 

X-« > • • • > U, > • • . • J U , ••••> X 

1 pi pi p 


can be written 


B. 


u 


1’ 


» ✓ J U X 

p p “ 

X ^ X . 

p P-1 


Af . T^ . . 

4 ^ ^1’ ' • • p ^1 


If we express both arrays as 1 dimensional vectors A and B , for which 

indexes are obtained by lexicographic ordering of the indexes (i^ , . . . . , i ) 

1 P 

and (u^, ..., multidimensional transfom can be expressed as a 

1-dimensional transform: 
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A = [[T^]@ [T^J ... ® [T^]] B 

A = [T] B 

The multidimensional transform has been reduced to a 1 dimensional trans- 
form. This expression now allows the evalution of the number of elementary 
operations and other generalizations discussed previously. 

7.3 Relations between transforms 

Two transforms with similar structures will often be related by matrix 
relations or energy invariants between the two sets of transformed 
coefficients . 

a) matrix relations between transforms of same order ; 

In [6] , matrix relations between the W-H and Haar transforms were 
proved. More generally, for WFH families, similar relations hold for all 
transforms lying on the same vertical line in the g-h graph of Fig. 9. 

These transforms only differ by the number of parent matrices [F 2 ] they 
include. Therefore a multiplication by all the missing [F^] matrices will 
generate one transform from the other. Note that these relations only 
involve computations in zones as defined in 3.3 or subzones (zonal 
divisions of a zone) . 

b) energy invariants ; 

By Parseval's theorem the total energy of the transform coefficients 
of a same vector with different transforms is preserved. However, it may 
happen that the energy of a subset of coefficients is the same for some trans 
forms: we say then there is an energy Invariant between these transforms. 

Energy invariants are most likely when the transforms have an identical struc 
ture with different multiplicators. For example, by direct comparison of the 






algorithms for the Fourier (Fig. 3b) W-H (Fig. 4c) and modified Haar 
(Fig. 5c), it is clear that the transformed coefficients before respective 
reorderings have identical energies in the zones defined in 3.3. This 
leads to the following energy invariants for the order 8. 


Zone 

Fourier 

(frequencies) 

W-H 

(zequencies) 

Mod . Haar 
(rank) 

0 

0 

0 

0 

1 

4 

7 

1 

2 

2,-2 

3,4 

2,3 

3 

1,3, -1,-3 

1,2, 5, 6 

4, 5,6, 7 


For the WFH families , the transforms with same sets of invariants 
form nested triangles as shown in Fig. 9: the introduction of additional 

factors leads to additional smaller subsets of coefficients of a same 
subzone over which energy is invariant; the relations between transforms which 
exist along vertical lines of the diagram of Fig. 9 preserves this energy 
invariance in zones. The invariants between the Generalized Discrete 
Transforms and the Modified Discrete Transforms have been studied by 
Rao et al. [18]. 

Conclusions 

In this work we have presented a unified treatment of unitary 
transforms having a fast algorithm. The use of recursive rules to describe 
unitary transforms allows a systematic way to view known transforms, 
to generate new transforms and provide a general approach to the evaluation 
of the computational complexity of transform algorithms. Among transforms 
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which are clearly related, we have studied the IC^ families and the WFH 
subfamilies which include most of the transforms considered in the 
literature. 

In addition to allowing the introduction of new transforms with 
properties of interest, the framework provided can be used in several 
other studies and applications of unitary transforms. In particular an 
error analysis of unitary transforms is being carried out. 
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